#Nature's Kidneys: the Role of Wetland Reserve Easements in Restoring Water Quality
#Nicole Karwowski and Marin Skidmore
#Replication File
#2025

#Supplementary Tables and Figures

#Packages
library(lubridate)
library(ggplot2)
library(dplyr)
library(fixest)
library(plm)
library(ggplot2)
library(readr)
library(mapview)
library(stringr)
library(stargazer)
library(cowplot)
library(broom)
library(openxlsx)
library(boot)
library(pvaluefunctions)
library(ggfixest)
library(patchwork)
library(Cairo)


#set wd
setwd("C:/Users/n57m645/OneDrive - Montana State University/UW_Madison/Water quality/Replication Package")
#open data

#load the final dataset 
load("Data/Main/main_huc12_panel.RData")


#keep only relevant years 
merge<-subset(merge, as.numeric(Year)>=1985)

#keep only observations that have an ammonia or phosphorous measurement
merge<-subset(merge, !is.na(ammonia_filtered) | !is.na(TP_filtered) | !is.na(TKN_filtered))

#give ever-treated their own df 
ever_treated <- subset(merge, merge$EverTreated==1)
set.seed(10)


#set variables corresponding to weather and agricultural quantiles 
breaks <- quantile(ever_treated$Ppt, probs = c(0.1, 0.33, 0.66, 0.9), na.rm = TRUE )
merge$PPT_bin <- cut(merge$Ppt, breaks = c(-Inf, breaks, Inf), 
                     labels = c("Very Low", "Low", "Medium", "High", "VeryHigh"), 
                     include.lowest = TRUE)
print(quantile(ever_treated$Ppt, probs = c(0.1, 0.33, 0.66, 0.9)))

breaks<- quantile(ever_treated$Temp, probs = c(0.1, 0.33, 0.66, 0.9), na.rm = TRUE )
merge$Temp_bin <- cut(merge$Temp, breaks = c(-Inf, breaks, Inf), 
                      labels = c("Very Low", "Low", "Medium", "High", "VeryHigh"), 
                      include.lowest = TRUE)
print(quantile(ever_treated$Temp, probs = c(0.1, 0.33, 0.66, 0.9)))
#################################################################################

#Create variable labels for tables
fixest::setFixest_dict(c(ammonia_filtered="Ammonia",
                         TP_filtered="Phosphorus",
                         TKN_filtered="Total Kjedahl Nitrogen",
                         RestoredBinary="Restored wetland (0/1)",
                         RestoredCumP="Restored wetland (prop)",
                         RestoredCum="Restored wetland (acres)",
                         TP_filtered_ihs="IHS(phosphorus)",
                         ammonia_filtered_ihs="IHS(ammonia)",
                         TKN_filtered_ihs="IHS(TKN)",
                         ihs_restored="IHS restored wetland (acres)",
                         ihs_upstream_restored_sf="IHS upstream restored wetland (acres)",
                         RestoredCumP_ihs="IHS restored wetland (prop)",
                         RestoredCum_ihs="IHS restored wetland (acres)",
                         UpstreamRestoredCumSfP="Upstream restored wetland (prop)",
                         ihs_upstream_restored_sf="IHS upstream restored wetland (acres)",
                         UpstreamRestoredCumSf="Upstream restored wetland (acres)",
                         RestoredBinaryUpstream="Upstream restored wetland (01/)",
                         UpstreamAmmoniaFiltered="Upstream ammonia",
                         UpstreamTpFiltered="Upstream phosphorus",
                         UpstreamTknFiltered="Upstream TKN",
                         upstream_sf_ammonia_filtered="Upstream ammonia (w=streamflow)",
                         upstream_sf_TP_filtered="Upstream phosphorus (w=streamflow)",
                         upstream_sf_TKN_filtered="Upstream TKN (w=streamflow)",
                         ihs_upstream_ammonia_sf = "Upstream ammonia (w=streamflow)",
                         ihs_upstream_TKN_sf = "Upstream TKN (w=streamflow)",
                         ihs_upstream_TP_sf = "Upstream phosphorus (w=streamflow)",
                         Ppt="Ppt",
                         PPT_bin="Ppt",
                         Temp="Temp",
                         Temp_bin="Temp",
                         barren="Barren",
                         cropland="Cropland",
                         developed="Developed",
                         vegetation="Vegetation",
                         water="Water",
                         grass_shrub="Grass shrub",
                         ice_snow="Ice snow",
                         tree_cover="Tree cover",
                         AnimalUnits="Animal units",
                         tri_ammonia="Ammonia dischage",
                         tri_ammonia_water="Ammonia discharge (water)",
                         tri_nitrate="Nitrate discharge",
                         tri_nitrate_water="Nitrate discharge (water)", 
                         tri_phosphorus="Phosphorus discharge",
                         tri_phosphorus_water="Phosphorus discharge (water)",
                         huc8="HUC8",
                         huc4="HUC4",
                         huc12="HUC12",
                         huc4year="HUC4 x year",
                         huc4month="HUC4 x month",
                         Year="Year",
                         Month="Month",
                         TTTyear2="Period",
                         TTT.closedyear2 = "Period",
                         EverTreated="",
                         times=""))
#fixed effect model
#set up macros
setFixest_fml(..ctrl_n = ~ +Ppt+Temp+cropland+developed+              
                vegetation+water + tri_ammonia_water,
              ..ctrl_p = ~ +Ppt+Temp+cropland+developed+              
                vegetation+water + tri_phosphorus_water,
              ..huc8_fe =~huc8 + Year + Month,
              ..huc12_fe = ~huc12 + Year + Month,
              ..huc4year_fe = ~huc8 + huc4year + Month,
              ..int_n = ~+PPT_bin + Temp_bin + cropland+developed+              
                vegetation+water + tri_ammonia_water +
                ihs_restored*PPT_bin+ihs_restored*Temp_bin,
              ..int_p = ~+PPT_bin + Temp_bin + cropland+developed+              
                vegetation+water + tri_phosphorus_water +
               ihs_restored*PPT_bin+ihs_restored*Temp_bin)


#post-processing functions for table formatting 
remove_model_numbers <- function(x) {
  # Remove lines starting specifically with & (1) & (2)
  x <- gsub("^\\s*& \\(1\\) & \\(2\\).*\n", "", x, perl = TRUE)
  return(x)
}

change_column_widths<- function(x) {
  #adjust the column widths of the table
  x <- gsub("{lcccccc}", "{m{6cm} m{1.4cm} m{1.4cm} m{1.4cm} m{1.4cm} m{1.4cm} m{1.4cm}}", x, perl=TRUE)
  return(x)
}

#function to remove model numbers and adjust column width
both<-function(x){
  x <- gsub("^\\s*& \\(1\\) & \\(2\\).*\n", "", x, perl = TRUE)
  x <- gsub("{lcccccc}", "{m{6cm} m{1.4cm} m{1.4cm} m{1.4cm} m{1.4cm} m{1.4cm} m{1.4cm}}", x, perl=TRUE)
  return(x)
}

#function to count the number of huc12 in our models
count_huc12 <- function(model, data, unit_var = "huc12") {
  # Extract observation row indices from the model
  obs_model <- obs(model)  
  
  # Ensure the dataset has a row index
  if (!"row_id" %in% names(data)) {
    data$row_id <- 1:nrow(data)  
  }
  
  # Filter dataset to match rows used in the model
  filtered_data <- data[data$row_id %in% obs_model, ]
  
  # Count unique units in the specified column
  return(length(unique(filtered_data[[unit_var]])))
}

#Function to county the number of unique HUCs in the control group
count_control_huc12 <- function(model, data, ref_years, treat_var = "first.treat.year", unit_var = "huc12") {
  # Get estimation sample
  obs_model <- obs(model)
  
  # Add row ID if needed
  if (!"row_id" %in% names(data)) {
    data$row_id <- seq_len(nrow(data))
  }
  
  # Filter data to estimation sample
  filtered_data <- data[data$row_id %in% obs_model, ]
  
  # Subset to control units (i.e., units whose treatment year is in ref_years)
  control_units <- filtered_data[filtered_data[[treat_var]] %in% ref_years, ]
  
  # Return number of unique HUC12s in control group
  return(length(unique(control_units[[unit_var]])))
}



###############################################################################################
#TWFE models in order of paper 

#model labels start with a for ammonia, p for phosphorus, n for nitrogen 
#start with nutrient for level outcome, ihs for ihs outcome 
#1 includes huc8, year, month fe
#2 includes huc8, huc4 x year, month fe
#3 includes huc12, year, month fe 
#control for upstream wq 

a1u_wq<-feols(fml=ammonia_filtered~ihs_restored+upstream_sf_ammonia_filtered +..ctrl_n| ..huc8_fe,
              cluster="huc8year", data=merge, subset=merge$EverTreated==1)

a2u_wq<-feols(fml=ammonia_filtered~ihs_restored+upstream_sf_ammonia_filtered+..ctrl_n| ..huc4year_fe,
              cluster="huc8year", data=merge, subset=merge$EverTreated==1)

a3u_wq<-feols(fml=ammonia_filtered~ihs_restored+upstream_sf_ammonia_filtered+..ctrl_n| ..huc12_fe,
              cluster="huc8year", data=merge, subset=merge$EverTreated==1)

ihsa1u_wq<-feols(fml=ihs_ammonia~ihs_restored+ihs_upstream_ammonia_sf+..ctrl_n| ..huc8_fe,
                 cluster="huc8year", data=merge, subset=merge$EverTreated==1)

ihsa2u_wq<-feols(fml=ihs_ammonia~ihs_restored+ihs_upstream_ammonia_sf+..ctrl_n| ..huc4year_fe,
                 cluster="huc8year", data=merge, subset=merge$EverTreated==1)

ihsa3u_wq<-feols(fml=ihs_ammonia~ihs_restored+ihs_upstream_ammonia_sf+..ctrl_n| ..huc12_fe,
                 cluster="huc8year", data=merge, subset=merge$EverTreated==1)

p1u_wq<-feols(fml=TP_filtered~ihs_restored+upstream_sf_TP_filtered+..ctrl_p| ..huc8_fe,
              cluster="huc8year", data=merge, subset=merge$EverTreated==1)

p2u_wq<-feols(fml=TP_filtered~ihs_restored+upstream_sf_TP_filtered+..ctrl_p| ..huc4year_fe,
              cluster="huc8year", data=merge, subset=merge$EverTreated==1)

p3u_wq<-feols(fml=TP_filtered~ihs_restored+upstream_sf_TP_filtered+..ctrl_p| ..huc12_fe,
              cluster="huc8year", data=merge, subset=merge$EverTreated==1)

ihsp1u_wq<-feols(fml=ihs_TP~ihs_restored+ihs_upstream_TP_sf+..ctrl_p| ..huc8_fe,
                 cluster="huc8year", data=merge, subset=merge$EverTreated==1)

ihsp2u_wq<-feols(fml=ihs_TP~ihs_restored+ihs_upstream_TP_sf+..ctrl_p| ..huc4year_fe,
                 cluster="huc8year", data=merge, subset=merge$EverTreated==1)

ihsp3u_wq<-feols(fml=ihs_TP~ihs_restored+ihs_upstream_TP_sf+..ctrl_p| ..huc12_fe,
                 cluster="huc8year", data=merge, subset=merge$EverTreated==1)

n1u_wq<-feols(fml=TKN_filtered~ihs_restored+upstream_sf_TKN_filtered+..ctrl_n| ..huc8_fe,
              cluster="huc8year", data=merge, subset=merge$EverTreated==1)

n2u_wq<-feols(fml=TKN_filtered~ihs_restored+upstream_sf_TKN_filtered+..ctrl_n| ..huc4year_fe,
              cluster="huc8year", data=merge, subset=merge$EverTreated==1)

n3u_wq<-feols(fml=TKN_filtered~ihs_restored+upstream_sf_TKN_filtered+..ctrl_n| ..huc12_fe,
              cluster="huc8year", data=merge, subset=merge$EverTreated==1)

ihsn1u_wq<-feols(fml=ihs_TKN~ihs_restored+ihs_upstream_TKN_sf+..ctrl_n| ..huc8_fe,
                 cluster="huc8year", data=merge, subset=merge$EverTreated==1)

ihsn2u_wq<-feols(fml=ihs_TKN~ihs_restored+ihs_upstream_TKN_sf+..ctrl_n| ..huc4year_fe,
                 cluster="huc8year", data=merge, subset=merge$EverTreated==1)

ihsn3u_wq<-feols(fml=ihs_TKN~ihs_restored+ihs_upstream_TKN_sf+..ctrl_n| ..huc12_fe,
                 cluster="huc8year", data=merge, subset=merge$EverTreated==1)

#upstream wq table 
fixest::etable(a1u_wq, a2u_wq, a3u_wq, ihsa1u_wq, ihsa2u_wq, ihsa3u_wq, 
               depvar = FALSE, drop.section = "fixef",
               tex=TRUE, digits=3, fitstat=c("my", "n", "ar2"), 
               keep=c("Restored", "IHS", "Upstream"),
               file="Supplementary Appendix/ammonia_treatedhuc12_upstreamWQ.tex",
               replace = TRUE,
               postprocess.tex = change_column_widths,
               style.tex=fixest::style.tex("aer", line.bottom = "",
                                           var.title = "$\\textbf{Panel A: ammonia}$"),
               extralines=list("__ "=c(" ", " ", " ", " ", " ", " "),
                               "^^  "=c(" ", " ", " ", " ", " ", " ")))

fixest::etable(n1u_wq, n2u_wq, n3u_wq, ihsn1u_wq, ihsn2u_wq, ihsn3u_wq, 
               depvar = FALSE, drop.section = "fixef",
               tex=TRUE, digits=3, fitstat=c("my", "n", "ar2"), 
               keep=c("Restored", "IHS", "Upstream"),
               file="Supplementary Appendix/tkn_treatedhuc12_upstreamWQ.tex",
               replace = TRUE,
               postprocess.tex = both,
               style.tex=fixest::style.tex("aer", line.top = "", line.bottom = "",
                                           var.title = "$\\textbf{Panel B: TKN}$"),
               extralines=list("__ "=c(" ", " ", " ", " ", " ", " "),
                               "^^  "=c(" ", " ", " ", " ", " ", " ")))

fixest::etable(p1u_wq, p2u_wq, p3u_wq, ihsp1u_wq, ihsp2u_wq, ihsp3u_wq, 
               depvar = FALSE, drop.section = "fixef",
               tex=TRUE, digits=3, fitstat=c("my", "n", "ar2"), 
               keep=c("Restored", "IHS", "Upstream"),
               file="Supplementary Appendix/phosphorus_treatedhuc12_upstreamWQ.tex",
               replace = TRUE,
               style.tex=fixest::style.tex("aer", line.top ="",
                                           fixef.title = "\\midrule",
                                           var.title = "$\\textbf{Panel C: phosphorus}$"),
               postprocess.tex = both,
               extralines=list("__"=c(" ", " ", " ", " "," ", " "),
                               "__Ever treated"=c("$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$","$\\checkmark$", "$\\checkmark$"),
                               "__Controls"=c("$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$","$\\checkmark$", "$\\checkmark$"),
                               "__"=c(" ", " ", " ", " "," ", " "),
                               "__HUC8 fixed effects"=c("$\\checkmark$",  "$\\checkmark$"," ", "$\\checkmark$", "$\\checkmark$"," "),
                               "__HUC12 fixed effects"=c(" "," ", "$\\checkmark$", " ", " ", "$\\checkmark$"),
                               "__Year fixed effects"=c("$\\checkmark$"," ", "$\\checkmark$", "$\\checkmark$"," ","$\\checkmark$"),
                               "__HUC4 x year fixed effects"=c(" ", "$\\checkmark$", " "," ", "$\\checkmark$"," "),
                               "__Month fixed effects"=c("$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$","$\\checkmark$", "$\\checkmark$")))


#all HUC models 
a1A<-feols(fml=ammonia_filtered~ihs_restored+..ctrl_n | ..huc8_fe,
           cluster="huc8year", data=merge)

a2A<-feols(fml=ammonia_filtered~ihs_restored+..ctrl_n | ..huc4year_fe,
           cluster="huc8year", data=merge)

a3A<-feols(fml=ammonia_filtered~ihs_restored+..ctrl_n| ..huc12_fe,
           cluster="huc8year", data=merge)

a1uA<-feols(fml=ammonia_filtered~ihs_restored+ihs_upstream_restored_sf+..ctrl_n | ..huc8_fe,
            cluster="huc8year", data=merge)

a2uA<-feols(fml=ammonia_filtered~ihs_restored+ihs_upstream_restored_sf+..ctrl_n | ..huc4year_fe,
            cluster="huc8year", data=merge)

a3uA<-feols(fml=ammonia_filtered~ihs_restored+ihs_upstream_restored_sf+..ctrl_n | ..huc12_fe,
            cluster="huc8year", data=merge)

p1A<-feols(fml=TP_filtered~ihs_restored+..ctrl_p | ..huc8_fe,
           cluster="huc8year", data=merge)

p2A<-feols(fml=TP_filtered~ihs_restored+..ctrl_p| ..huc4year_fe,
           cluster="huc8year", data=merge)

p3A<-feols(fml=TP_filtered~ihs_restored+..ctrl_p | ..huc12_fe,
           cluster="huc8year", data=merge)

p1uA<-feols(fml=TP_filtered~ihs_restored+ihs_upstream_restored_sf+..ctrl_p | ..huc8_fe,
            cluster="huc8year", data=merge)

p2uA<-feols(fml=TP_filtered~ihs_restored+ihs_upstream_restored_sf+..ctrl_p | ..huc4year_fe,
            cluster="huc8year", data=merge)

p3uA<-feols(fml=TP_filtered~ihs_restored+ihs_upstream_restored_sf+..ctrl_p | ..huc12_fe,
            cluster="huc8year", data=merge)


n1A<-feols(fml=TKN_filtered~ihs_restored+..ctrl_n | ..huc8_fe,
           cluster="huc8year", data=merge)

n2A<-feols(fml=TKN_filtered~ihs_restored+..ctrl_n| ..huc4year_fe,
           cluster="huc8year", data=merge)

n3A<-feols(fml=TKN_filtered~ihs_restored+..ctrl_n | ..huc12_fe,
           cluster="huc8year", data=merge)

n1uA<-feols(fml=TKN_filtered~ihs_restored+ihs_upstream_restored_sf+..ctrl_n | ..huc8_fe,
            cluster="huc8year", data=merge)

n2uA<-feols(fml=TKN_filtered~ihs_restored+ihs_upstream_restored_sf+..ctrl_n | ..huc4year_fe,
            cluster="huc8year", data=merge)

n3uA<-feols(fml=TKN_filtered~ihs_restored+ihs_upstream_restored_sf+..ctrl_n | ..huc12_fe,
            cluster="huc8year", data=merge)


#all HUC table 
fixest::etable(a1A, a2A, a3A, a1uA, a2uA, a3uA, 
               depvar = FALSE, drop.section="fixef",
               tex=TRUE, digits=3, fitstat=c("my", "n", "ar2"), 
               keep=c("IHS"),
               file="Supplementary Appendix/ammonia_allhuc12_FE.tex",
               replace = TRUE,               
               postprocess.tex = change_column_widths,
               style.tex=fixest::style.tex("aer", line.bottom = "",
                                           var.title = "$\\textbf{Panel A: ammonia}$"),
               extralines=list("__ "=c(" ", " ", " ", " ", " ", " "),
                               "^^  "=c(" ", " ", " ", " ", " ", " ")))
fixest::etable(n1A, n2A, n3A, n1uA, n2uA, n3uA, 
               depvar = FALSE, drop.section="fixef",
               tex=TRUE, digits=3, fitstat=c("my", "n", "ar2"), 
               keep=c("IHS"),
               file="Supplementary Appendix/tkn_allhuc12_FE.tex",
               replace = TRUE,               
               postprocess.tex = both,
               style.tex=fixest::style.tex("aer", line.top = "", line.bottom = "",
                                           var.title = "$\\textbf{Panel B: TKN}$"),
               extralines=list("__ "=c(" ", " ", " ", " ", " ", " "),
                               "^^  "=c(" ", " ", " ", " ", " ", " ")))


fixest::etable(p1A, p2A, p3A, p1uA, p2uA, p3uA,
               depvar = FALSE, drop.section = "fixef",
               tex=TRUE, digits=3, fitstat=c("my", "n", "ar2"), 
               keep=c("IHS"),
               file="Supplementary Appendix/phosphorus_allhuc12_FE.tex",
               replace = TRUE,
               style.tex=fixest::style.tex("aer", line.top ="",
                                           fixef.title = "\\midrule",
                                           var.title = "$\\textbf{Panel C: phosphorus}$"),
               postprocess.tex = both,                
               extralines=list("__"=c(" ", " ", " ", " "," ", " "),
                               "__Ever treated"=c("$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$","$\\checkmark$", "$\\checkmark$"),
                               "__Controls"=c("$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$","$\\checkmark$", "$\\checkmark$"),
                               "__"=c(" ", " ", " ", " "," ", " "),
                               "__HUC8 fixed effects"=c("$\\checkmark$",  "$\\checkmark$"," ", "$\\checkmark$", "$\\checkmark$"," "),
                               "__HUC12 fixed effects"=c(" "," ", "$\\checkmark$", " ", " ", "$\\checkmark$"),
                               "__Year fixed effects"=c("$\\checkmark$"," ", "$\\checkmark$", "$\\checkmark$"," ","$\\checkmark$"),
                               "__HUC4 x year fixed effects"=c(" ", "$\\checkmark$", " "," ", "$\\checkmark$"," "),
                               "__Month fixed effects"=c("$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$","$\\checkmark$", "$\\checkmark$")))

#IHS WQ outcome of main models 

ihsa1<-feols(fml=ihs_ammonia~ihs_restored+..ctrl_n | ..huc8_fe, 
             cluster="huc8year", data=merge, subset=merge$EverTreated==1)

ihsa2<-feols(fml=ihs_ammonia~ihs_restored+..ctrl_n| ..huc4year_fe,
             cluster="huc8year", data=merge, subset=merge$EverTreated==1)

ihsa3<-feols(fml=ihs_ammonia~ihs_restored+..ctrl_n| ..huc12_fe,
             cluster="huc8year", data=merge, subset=merge$EverTreated==1)

ihsp1<-feols(fml=ihs_TP~ihs_restored+..ctrl_p | ..huc8_fe, 
             cluster="huc8year", data=merge, subset=merge$EverTreated==1)

ihsp2<-feols(fml=ihs_TP~ihs_restored+..ctrl_p| ..huc4year_fe,
             cluster="huc8year", data=merge, subset=merge$EverTreated==1)

ihsp3<-feols(fml=ihs_TP~ihs_restored+..ctrl_p| ..huc12_fe,
             cluster="huc8year", data=merge, subset=merge$EverTreated==1)

ihsn1<-feols(fml=ihs_TKN~ihs_restored+..ctrl_n | ..huc8_fe, 
             cluster="huc8year", data=merge, subset=merge$EverTreated==1)

ihsn2<-feols(fml=ihs_TKN~ihs_restored+..ctrl_n| ..huc4year_fe,
             cluster="huc8year", data=merge, subset=merge$EverTreated==1)

ihsn3<-feols(fml=ihs_TKN~ihs_restored+..ctrl_n| ..huc12_fe,
             cluster="huc8year", data=merge, subset=merge$EverTreated==1)



#ihs WQ plus upstream restored 
ihsa1u<-feols(fml=ihs_ammonia~ihs_restored+ihs_upstream_restored_sf+..ctrl_n| ..huc8_fe,
              cluster="huc8year", data=merge, subset=merge$EverTreated==1)

ihsa2u<-feols(fml=ihs_ammonia~ihs_restored+ihs_upstream_restored_sf+..ctrl_n| ..huc4year_fe,
              cluster="huc8year", data=merge, subset=merge$EverTreated==1)

ihsa3u<-feols(fml=ihs_ammonia~ihs_restored+ihs_upstream_restored_sf+..ctrl_n| ..huc12_fe,
              cluster="huc8year", data=merge, subset=merge$EverTreated==1)

ihsp1u<-feols(fml=ihs_TP~ihs_restored+ihs_upstream_restored_sf+..ctrl_p| ..huc8_fe,
              cluster="huc8year", data=merge, subset=merge$EverTreated==1)

ihsp2u<-feols(fml=ihs_TP~ihs_restored+ihs_upstream_restored_sf+..ctrl_p| ..huc4year_fe,
              cluster="huc8year", data=merge, subset=merge$EverTreated==1)

ihsp3u<-feols(fml=ihs_TP~ihs_restored+ihs_upstream_restored_sf+..ctrl_p| ..huc12_fe,
              cluster="huc8year", data=merge, subset=merge$EverTreated==1)

ihsn1u<-feols(fml=ihs_TKN~ihs_restored+ihs_upstream_restored_sf+..ctrl_n| ..huc8_fe,
              cluster="huc8year", data=merge, subset=merge$EverTreated==1)

ihsn2u<-feols(fml=ihs_TKN~ihs_restored+ihs_upstream_restored_sf+..ctrl_n| ..huc4year_fe,
              cluster="huc8year", data=merge, subset=merge$EverTreated==1)

ihsn3u<-feols(fml=ihs_TKN~ihs_restored+ihs_upstream_restored_sf+..ctrl_n| ..huc12_fe,
              cluster="huc8year", data=merge, subset=merge$EverTreated==1)

#IHS WQ table 
fixest::etable(ihsa1, ihsa2, ihsa3, ihsa1u, ihsa2u, ihsa3u, 
               depvar = FALSE, drop.section = "fixef",
               tex=TRUE, digits=3, fitstat=c("my", "n", "ar2"), 
               keep=c("IHS", "IHS"),
               file="Supplementary Appendix/ihs_ammonia_treatedhuc12_FE.tex",
               replace = TRUE,
               postprocess.tex = change_column_widths,
               style.tex=fixest::style.tex("aer", line.bottom = "",
                                           var.title = "$\\textbf{Panel A: ammonia}$"),
               extralines=list("__ "=c(" ", " ", " ", " ", " ", " "),
                               "^^  "=c(" ", " ", " ", " ", " ", " ")))

fixest::etable(ihsn1, ihsn2, ihsn3, ihsn1u, ihsn2u, ihsn3u, 
               depvar = FALSE, drop.section = "fixef",
               tex=TRUE, digits=3, fitstat=c("my", "n", "ar2"), 
               keep=c("IHS", "IHS"),
               file="Supplementary Appendix/ihs_tkn_treatedhuc12_FE.tex",
               replace = TRUE,
               postprocess.tex = both,
               style.tex=fixest::style.tex("aer", line.top = "", line.bottom = "",
                                           var.title = "$\\textbf{Panel B: TKN}$"),
               extralines=list("__ "=c(" ", " ", " ", " ", " ", " "),
                               "^^  "=c(" ", " ", " ", " ", " ", " ")))


fixest::etable(ihsp1, ihsp2, ihsp3, ihsp1u, ihsp2u, ihsp3u, 
               depvar = FALSE, drop.section = "fixef",
               tex=TRUE, digits=3, fitstat=c("my", "n", "ar2"), 
               keep=c("IHS"),
               file="Supplementary Appendix/ihs_phosphorus_treatedhuc12_FE.tex",
               replace = TRUE,
               postprocess.tex = both,
               style.tex=fixest::style.tex("aer", line.top ="",
                                           fixef.title = "\\midrule",
                                           var.title = "$\\textbf{Panel C: phosphorus}$"),
               extralines=list("__"=c(" ", " ", " ", " "," ", " "),
                               "__Ever treated"=c("$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$","$\\checkmark$", "$\\checkmark$"),
                               "__Controls"=c("$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$","$\\checkmark$", "$\\checkmark$"),
                               "__"=c(" ", " ", " ", " "," ", " "),
                               "__HUC8 fixed effects"=c("$\\checkmark$",  "$\\checkmark$"," ", "$\\checkmark$", "$\\checkmark$"," "),
                               "__HUC12 fixed effects"=c(" "," ", "$\\checkmark$", " ", " ", "$\\checkmark$"),
                               "__Year fixed effects"=c("$\\checkmark$"," ", "$\\checkmark$", "$\\checkmark$"," ","$\\checkmark$"),
                               "__HUC4 x year fixed effects"=c(" ", "$\\checkmark$", " "," ", "$\\checkmark$"," "),
                               "__Month fixed effects"=c("$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$","$\\checkmark$", "$\\checkmark$")))



#percent wetland 

a1_p<-feols(fml=ammonia_filtered~RestoredCumP+..ctrl_n | ..huc8_fe, 
            cluster="huc8year", data=merge, subset=merge$EverTreated==1)

a2_p<-feols(fml=ammonia_filtered~RestoredCumP+..ctrl_n| ..huc4year_fe,
            cluster="huc8year", data=merge, subset=merge$EverTreated==1)

a3_p<-feols(fml=ammonia_filtered~RestoredCumP+..ctrl_n| ..huc12_fe,
            cluster="huc8year", data=merge, subset=merge$EverTreated==1)

p1_p<-feols(fml=TP_filtered~RestoredCumP+..ctrl_p | ..huc8_fe, 
            cluster="huc8year", data=merge, subset=merge$EverTreated==1)

p2_p<-feols(fml=TP_filtered~RestoredCumP+..ctrl_p| ..huc4year_fe,
            cluster="huc8year", data=merge, subset=merge$EverTreated==1)

p3_p<-feols(fml=TP_filtered~RestoredCumP+..ctrl_p| ..huc12_fe,
            cluster="huc8year", data=merge, subset=merge$EverTreated==1)

n1_p<-feols(fml=TKN_filtered~RestoredCumP+..ctrl_n | ..huc8_fe, 
            cluster="huc8year", data=merge, subset=merge$EverTreated==1)

n2_p<-feols(fml=TKN_filtered~RestoredCumP+..ctrl_n| ..huc4year_fe,
            cluster="huc8year", data=merge, subset=merge$EverTreated==1)

n3_p<-feols(fml=TKN_filtered~RestoredCumP+..ctrl_n| ..huc12_fe,
            cluster="huc8year", data=merge, subset=merge$EverTreated==1)


#binary wetland 

a1_bin <-feols(fml=ammonia_filtered~RestoredBinary+..ctrl_n | ..huc8_fe, 
               cluster="huc8year", data=merge, subset=merge$EverTreated==1)

a2_bin<-feols(fml=ammonia_filtered~RestoredBinary+..ctrl_n| ..huc4year_fe,
              cluster="huc8year", data=merge, subset=merge$EverTreated==1)

a3_bin<-feols(fml=ammonia_filtered~RestoredBinary+..ctrl_n| ..huc12_fe,
              cluster="huc8year", data=merge, subset=merge$EverTreated==1)

p1_bin <-feols(fml=TP_filtered~RestoredBinary+..ctrl_p | ..huc8_fe, 
               cluster="huc8year", data=merge, subset=merge$EverTreated==1)

p2_bin<-feols(fml=TP_filtered~RestoredBinary+..ctrl_p| ..huc4year_fe,
              cluster="huc8year", data=merge, subset=merge$EverTreated==1)

p3_bin<-feols(fml=TP_filtered~RestoredBinary+..ctrl_p| ..huc12_fe,
              cluster="huc8year", data=merge, subset=merge$EverTreated==1)

n1_bin <-feols(fml=TKN_filtered~RestoredBinary+..ctrl_n | ..huc8_fe, 
               cluster="huc8year", data=merge, subset=merge$EverTreated==1)

n2_bin<-feols(fml=TKN_filtered~RestoredBinary+..ctrl_n| ..huc4year_fe,
              cluster="huc8year", data=merge, subset=merge$EverTreated==1)

n3_bin<-feols(fml=TKN_filtered~RestoredBinary+..ctrl_n| ..huc12_fe,
              cluster="huc8year", data=merge, subset=merge$EverTreated==1)

fixest::etable(a1_p, a2_p, a3_p, a1_bin, a2_bin, a3_bin, 
               depvar = FALSE, drop.section="fixef",
               tex=TRUE, digits=3, fitstat=c("my", "n", "ar2"), 
               keep=c("Restored"),
               file="Supplementary Appendix/ammonia_treatedhuc12_pct.tex",
               replace = TRUE,               
               postprocess.tex = change_column_widths,
               style.tex=fixest::style.tex("aer", line.bottom = "",
                                           var.title = "$\\textbf{Panel A: ammonia}$"),
               extralines=list("__ "=c(" ", " ", " ", " ", " ", " "),
                               "^^  "=c(" ", " ", " ", " ", " ", " ")))
fixest::etable(n1_p, n2_p, n3_p, n1_bin, n2_bin, n3_bin, 
               depvar = FALSE, drop.section="fixef",
               tex=TRUE, digits=3, fitstat=c("my", "n", "ar2"), 
               keep=c("Restored"),
               file="Supplementary Appendix/tkn_treatedhuc12_pct.tex",
               replace = TRUE,               
               postprocess.tex = both,
               style.tex=fixest::style.tex("aer", line.top = "", line.bottom = "",
                                           var.title = "$\\textbf{Panel B: TKN}$"),
               extralines=list("__ "=c(" ", " ", " ", " ", " ", " "),
                               "^^  "=c(" ", " ", " ", " ", " ", " ")))


fixest::etable(p1_p, p2_p, p3_p, p1_bin, p2_bin, p3_bin,
               depvar = FALSE, drop.section = "fixef",
               tex=TRUE, digits=3, fitstat=c("my", "n", "ar2"), 
               keep=c("Restored"),
               file="Supplementary Appendix/phosphorus_treatedhuc12_pct.tex",
               replace = TRUE,
               style.tex=fixest::style.tex("aer", line.top ="",
                                           fixef.title = "\\midrule",
                                           var.title = "$\\textbf{Panel C: phosphorus}$"),
               postprocess.tex = both,
               extralines=list("__"=c(" ", " ", " ", " "," ", " "),
                               "__Ever treated"=c("$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$","$\\checkmark$", "$\\checkmark$"),
                               "__Controls"=c("$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$","$\\checkmark$", "$\\checkmark$"),
                               "__"=c(" ", " ", " ", " "," ", " "),
                               "__HUC8 fixed effects"=c("$\\checkmark$",  "$\\checkmark$"," ", "$\\checkmark$", "$\\checkmark$"," "),
                               "__HUC12 fixed effects"=c(" "," ", "$\\checkmark$", " ", " ", "$\\checkmark$"),
                               "__Year fixed effects"=c("$\\checkmark$"," ", "$\\checkmark$", "$\\checkmark$"," ","$\\checkmark$"),
                               "__HUC4 x year fixed effects"=c(" ", "$\\checkmark$", " "," ", "$\\checkmark$"," "),
                               "__Month fixed effects"=c("$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$","$\\checkmark$", "$\\checkmark$")))



##############################################################################################
#Event studies 
#appendix only has ever-treated HUC12 models 

#Ammonia 
#basic event study (i.e., no controls)
#ever treated
#TWFE
a1es_uc <-fixest::feols(ammonia_filtered~i(TTTyear2, EverTreated, ref=-1) | ..huc8_fe , 
                        data=merge, cluster="huc8year", subset=merge$EverTreated==1)

#sun and abraham 
a2es_uc<-fixest::feols(ammonia_filtered~sunab(first.treat.year, TTTyear2, ref.c = 2018) | ..huc8_fe, 
                       data=merge, cluster="huc8year", subset=merge$EverTreated==1 & merge$Year < 2018)

#conditional event study 
#TWFE
a1es<-fixest::feols(ammonia_filtered~i(TTTyear2, EverTreated, ref=-1) + ..ctrl_n | ..huc8_fe, 
                    data=merge, cluster="huc8year", subset=merge$EverTreated==1)

#sun and abraham 
a2es<-fixest::feols(ammonia_filtered~sunab(first.treat.year, TTTyear2, ref.c = 2018) + ..ctrl_n | ..huc8_fe, 
                    data=merge, cluster="huc8year", subset=merge$EverTreated==1& merge$Year < 2018)

#HUC12 FE 
#sun and abraham only
#with and without controls 
#ammonia
a2es_12fe<-fixest::feols(ammonia_filtered~sunab(first.treat.year, TTTyear2, ref.c = 2018) +..ctrl_n | ..huc12_fe, 
                         data=merge, cluster="huc8year", subset=merge$EverTreated==1& merge$Year < 2018)

#without controls 
a2es_uc_12fe<-fixest::feols(ammonia_filtered~sunab(first.treat.year, TTTyear2, ref.c = 2018)  | ..huc12_fe, 
                            data=merge, cluster="huc8year", subset=merge$EverTreated==1& merge$Year < 2018)

#sun and abraham using ihs water quality 
a2es_uc_ihs <-fixest::feols(ihs_ammonia ~sunab(first.treat.year, TTTyear2, ref.c = 2018) | ..huc8_fe, 
                            data=merge, cluster="huc8year", subset=merge$EverTreated==1& merge$Year < 2018)

a2es_ihs <-fixest::feols(ihs_ammonia ~sunab(first.treat.year, TTTyear2, ref.c = 2018) +..ctrl_n | ..huc8_fe, 
                         data=merge, cluster="huc8year", subset=merge$EverTreated==1& merge$Year < 2018)

#Contract closure 
#ammonia
#without controls
#treated
a2es_uc_closed<-fixest::feols(ammonia_filtered~sunab(first.treat.year, TTT.closedyear2, ref.c = 2018)| ..huc8_fe, 
                              data=merge, cluster="huc8year", subset=merge$EverTreated==1& merge$Year < 2018)
#with controls
#treated
a2es_closed<-fixest::feols(ammonia_filtered~sunab(first.treat.year, TTT.closedyear2, ref.c = 2018) +..ctrl_n| ..huc8_fe, 
                           data=merge, cluster="huc8year", subset=merge$EverTreated==1& merge$Year < 2018)


#ATT 
#with controls 
a2es_att<-fixest::feols(ammonia_filtered~sunab(first.treat.year, TTTyear2, att = TRUE, ref.c = 2018) +..ctrl_n | ..huc8_fe, 
                        data=merge, cluster="huc8year", subset=merge$EverTreated==1& merge$Year < 2018)
#without controls 
#ever treated
a2es_uc_att<-fixest::feols(ammonia_filtered~sunab(first.treat.year, TTTyear2, att = TRUE, ref.c = 2018)  | ..huc8_fe, 
                           data=merge, cluster="huc8year", subset=merge$EverTreated==1& merge$Year < 2018)


#Phosphorus
#basic event study
#ever treated
#TWFE
p1es_uc<-fixest::feols(TP_filtered~i(TTTyear2, EverTreated, ref=-1) | ..huc8_fe, 
                       data=merge, cluster="huc8year", subset=merge$EverTreated==1)

#sun and abraham 
p2es_uc<-fixest::feols(TP_filtered~sunab(first.treat.year, TTTyear2, ref.c = 2018) | ..huc8_fe, 
                       data=merge, cluster="huc8year", subset=merge$EverTreated==1& merge$Year < 2018)

#conditional event study -- main text
#TWFE
p1es<-fixest::feols(TP_filtered~i(TTTyear2, EverTreated, ref=-1)+..ctrl_p | ..huc8_fe, 
                    data=merge, cluster="huc8year", subset=merge$EverTreated==1)

#sun and abraham 
p2es<-fixest::feols(TP_filtered~sunab(first.treat.year, TTTyear2, ref.c = 2018)+..ctrl_p  | ..huc8_fe, 
                    data=merge, cluster="huc8year", subset=merge$EverTreated==1& merge$Year < 2018)

#HUC12 FE, SA only, with and without controls 
p2es_12fe<-fixest::feols(TP_filtered~sunab(first.treat.year, TTTyear2, ref.c = 2018)+..ctrl_p  | ..huc12_fe, 
                         data=merge, cluster="huc8year", subset=merge$EverTreated==1& merge$Year < 2018)

p2es_uc_12fe<-fixest::feols(TP_filtered~sunab(first.treat.year, TTTyear2, ref.c = 2018)  | ..huc12_fe, 
                            data=merge, cluster="huc8year", subset=merge$EverTreated==1& merge$Year < 2018)

#IHS WQ 

p2es_uc_ihs <-fixest::feols(ihs_TP ~sunab(first.treat.year, TTTyear2, ref.c = 2018) | ..huc8_fe, 
                            data=merge, cluster="huc8year", subset=merge$EverTreated==1& merge$Year < 2018)

p2es_ihs <-fixest::feols(ihs_TP ~sunab(first.treat.year, TTTyear2, ref.c = 2018) +..ctrl_p | ..huc8_fe, 
                         data=merge, cluster="huc8year", subset=merge$EverTreated==1& merge$Year < 2018)

#contract closure 
#without controls
#ever treated 
p2es_uc_closed<-fixest::feols(TP_filtered~sunab(first.treat.year, TTT.closedyear2, ref.c = 2018)  | ..huc8_fe, 
                              data=merge, cluster="huc8year", subset=merge$EverTreated==1& merge$Year < 2018)
#with controls
#ever treated
p2es_closed<-fixest::feols(TP_filtered~sunab(first.treat.year, TTT.closedyear2, ref.c = 2018)+..ctrl_p | ..huc8_fe, 
                           data=merge, cluster="huc8year", subset=merge$EverTreated==1& merge$Year < 2018)

#ATT
#with controls 
p2es_att<-fixest::feols(TP_filtered~sunab(first.treat.year, TTTyear2, att = TRUE, ref.c = 2018)+..ctrl_p  | ..huc8_fe, 
                        data=merge, cluster="huc8year", subset=merge$EverTreated==1& merge$Year < 2018)

#without controls
p2es_uc_att<-fixest::feols(TP_filtered~sunab(first.treat.year, TTTyear2, att = TRUE, ref.c = 2018)  | ..huc8_fe, 
                           data=merge, cluster="huc8year", subset=merge$EverTreated==1& merge$Year < 2018)



#TKN
#basic event study
#ever treated
#TWFE
n1es_uc<-fixest::feols(TKN_filtered~i(TTTyear2, EverTreated, ref=-1) | ..huc8_fe, 
                       data=merge, cluster="huc8year", subset=merge$EverTreated==1)

#sun and abraham 
n2es_uc<-fixest::feols(TKN_filtered~sunab(first.treat.year, TTTyear2, ref.c = 2018) | ..huc8_fe, 
                       data=merge, cluster="huc8year", subset=merge$EverTreated==1& merge$Year < 2018)

#conditional event study -- main text
#TWFE
n1es<-fixest::feols(TKN_filtered~i(TTTyear2, EverTreated, ref=-1)+..ctrl_n | ..huc8_fe, 
                    data=merge, cluster="huc8year", subset=merge$EverTreated==1)

#sun and abraham 
n2es<-fixest::feols(TKN_filtered~sunab(first.treat.year, TTTyear2, ref.c = 2018)+..ctrl_n  | ..huc8_fe, 
                    data=merge, cluster="huc8year", subset=merge$EverTreated==1& merge$Year < 2018)

#HUC12 FE, SA only, with and without controls 
n2es_12fe<-fixest::feols(TKN_filtered~sunab(first.treat.year, TTTyear2, ref.c = 2018)+..ctrl_n  | ..huc12_fe, 
                         data=merge, cluster="huc8year", subset=merge$EverTreated==1& merge$Year < 2018)

n2es_uc_12fe<-fixest::feols(TKN_filtered~sunab(first.treat.year, TTTyear2, ref.c = 2018)  | ..huc12_fe, 
                            data=merge, cluster="huc8year", subset=merge$EverTreated==1& merge$Year < 2018)

#IHS WQ 

n2es_uc_ihs <-fixest::feols(ihs_TKN ~sunab(first.treat.year, TTTyear2, ref.c = 2018) | ..huc8_fe, 
                            data=merge, cluster="huc8year", subset=merge$EverTreated==1& merge$Year < 2018)

n2es_ihs <-fixest::feols(ihs_TKN ~sunab(first.treat.year, TTTyear2, ref.c = 2018) +..ctrl_n | ..huc8_fe, 
                         data=merge, cluster="huc8year", subset=merge$EverTreated==1& merge$Year < 2018)

#contract closure 
#without controls
#ever treated 
n2es_uc_closed<-fixest::feols(TKN_filtered~sunab(first.treat.year, TTT.closedyear2, ref.c = 2018)  | ..huc8_fe, 
                              data=merge, cluster="huc8year", subset=merge$EverTreated==1& merge$Year < 2018)
#with controls
#ever treated
n2es_closed<-fixest::feols(TKN_filtered~sunab(first.treat.year, TTT.closedyear2, ref.c = 2018)+..ctrl_n | ..huc8_fe, 
                           data=merge, cluster="huc8year", subset=merge$EverTreated==1& merge$Year < 2018)


#ATT
#with controls 
n2es_att<-fixest::feols(TKN_filtered~sunab(first.treat.year, TTTyear2, att = TRUE, ref.c = 2018)+..ctrl_n  | ..huc8_fe, 
                        data=merge, cluster="huc8year", subset=merge$EverTreated==1& merge$Year < 2018)

#without controls
n2es_uc_att<-fixest::feols(TKN_filtered~sunab(first.treat.year, TTTyear2, att = TRUE, ref.c = 2018)  | ..huc8_fe, 
                           data=merge, cluster="huc8year", subset=merge$EverTreated==1 & merge$Year < 2018)



#Event study supplement graphs
#we're going to output these based on the format of the treatment variable 
ols_list <- list(
                "a1es" = "Ammonia\nOLS",
                "p1es" = "Phosphorus\nOLS",
                "a1es_uc" = "Ammonia\nOLS",
                 "p1es_uc" = "Phosphorus\nOLS",
                "n1es" = "Total Kjeldahl Nitrogen\nOLS",
                  "n1es_uc" = "Total Kjeldahl Nitrogen\nOLS")
  
sunab_list <- list(
                "a2es_uc" = "Ammonia\nS&A",
                "p2es_uc" = "Phosphorus\nS&A",
                "a2es_12fe" = "Ammonia\nControls",
                "a2es_uc_12fe" = "Ammonia\nNo controls",
                "p2es_12fe" = "Phosphorus\nControls",
                "p2es_uc_12fe" = "Phosphorus\nNo controls",
                "a2es_ihs" = "Ammonia\nControls",
                "a2es_uc_ihs" = "Ammonia\nNo controls",
                "p2es_ihs" = "Phosphorus\nControls",
                "p2es_uc_ihs" = "Phosphorus\nNo controls",
                "n2es_uc" = "Total Kjeldahl Nitrogen\nS&A",
                "n2es_12fe" = "Total Kjeldahl Nitrogen\nControls",
                "n2es_uc_12fe" = "Total Kjeldahl Nitrogen\nNo controls",
                "n2es_ihs" = "Total Kjeldahl Nitrogen\nControls",
                "n2es_uc_ihs" = "Total Kjeldahl Nitrogen\nNo controls")
  
  closed_list <- list(
                "a2es_closed" = "Ammonia\nControls",
                 "a2es_uc_closed" = "Ammonia\nNo controls",
                 "n2es_uc_closed" = "Total Kjeldahl Nitrogen\nNo controls",
                 "p2es_closed" = "Phosphorus\nControls",
                 "p2es_uc_closed" = "Phosphorus\nNo controls",
                 "n2es_closed" = "Total Kjeldahl Nitrogen\nControls")

model_groups <- list(
  A_and_P = c("a1es" ,
              "p1es",
              "a1es_uc",
              "p1es_uc",
              "a2es_uc" ,
              "p2es_uc" ,
              "a2es_12fe",
              "a2es_uc_12fe" ,
              "p2es_12fe" ,
              "p2es_uc_12fe" ,
              "a2es_ihs" ,
              "a2es_uc_ihs",
              "p2es_ihs" ,
              "p2es_uc_ihs",
              "a2es_closed",
              "a2es_uc_closed",
              "p2es_closed",
              "p2es_uc_closed"),
  TKN = c("n1es",
          "n1es_uc",
          "n2es_uc",
          "n2es_12fe",
          "n2es_uc_12fe",
          "n2es_ihs" ,
          "n2es_uc_ihs",
          "n2es_closed",
          "n2es_uc_closed")
)
#we also need conditional y axis breaks bc TKN has different distribution 
group_y_axis <- list(
  A_and_P = list(breaks = seq(-0.2, 0.2, by = 0.1), limits = c(-0.2, 0.2)),
  TKN = list(breaks = seq(-0.3, 0.3, by = 0.1), limits = c(-0.3, 0.3))
)
# use ggplot to make the event study graphs
My_Theme <- theme(
  title = element_text(size = 14),
  text = element_text(color = "black"),
  strip.text = element_text(color = "black", size = 14),
  strip.background = element_blank(),
  axis.title.x = element_text(size = 14),
  axis.text.x = element_text(size = 14),
  axis.text.y = element_text(size = 14),
  axis.title.y = element_text(size = 14),
  legend.position = "none",
  panel.background = element_rect(fill = "white", color = "black"),  # keeps white background with black border
  panel.grid.major = element_line(color = "gray85"),  # soft main grid lines
  panel.grid.minor = element_line(color = "gray90"),  # even softer minor lines
  plot.background = element_rect(fill = "white", color = NA)
)



#variable names
# var_labels<-as_labeller(c('ammonia_filtered'='Ammonia', 'TP_filtered'='Phosphorus', 'TKN_filtered'='Total kjedahl nitrogen',
#                           'All HUC12'= 'All HUC12', 'Ever-treated HUC12'= 'Ever-treated HUC12'))




for (model_name in names(ols_list)) {
  model <- get(model_name)  # retrieve the model object
  tidy_model <- broom::tidy(model)  # extract coefficients
  
  event_df <- tidy_model %>%
    filter(grepl("TTTyear2::.*:EverTreated", term)) %>%
      mutate(
        time = as.numeric(gsub("TTTyear2::(-?\\d+):?.*", "\\1", term)),
    ci_lower = estimate - 1.96 * std.error,
      ci_upper = estimate + 1.96 * std.error
    )
  event_df <- bind_rows(
    event_df,
    tibble(time = -1, estimate = 0, ci_lower = 0, ci_upper = 0)
  )

  # Identify group and pull settings
  current_group <- names(Filter(function(x) model_name %in% x, model_groups))[1]
  y_settings <- group_y_axis[[current_group]]
  
  # Plot
  # p <- ggplot(event_df, aes(x = time, y = estimate)) +
  #   geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper), fill = "gray80") +
  #   geom_line(color = "black") +
  #   geom_point(color = "black") +
  #   geom_vline(xintercept = -1, linetype = "dashed", color = "black") +
  #   labs(
  #     title = ols_list[[model_name]],
  #     x = "Event Time",
  #     y = "Estimate and 95% CI"
  #   ) +
  #   My_Theme+ 
  #   scale_y_continuous(breaks = y_settings$breaks,
  #                      labels = scales::label_number(accuracy = 0.01)) +
  #   coord_cartesian(ylim = y_settings$limits)  
  
  p <- ggplot(event_df, aes(x = time, y = estimate)) +
    geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper), fill = "gray60") +
    geom_line(color = "black") +
    geom_point(color = "black") +
    geom_vline(xintercept = -1, linetype = "dashed", color = "black") +
    geom_hline(yintercept = 0, linetype = "solid", color = "black") + 
    labs(
      title = ols_list[[model_name]],
      x = "Event Time",
      y = "Estimate and 95% CI"
    ) +
    My_Theme+ 
    scale_y_continuous(breaks = y_settings$breaks,
                       labels = scales::label_number(accuracy = 0.01)) +
    coord_cartesian(ylim = y_settings$limits)  +
    theme(plot.title = element_text(hjust = 0.5))
  
  
  filename <- paste0("Supplementary Appendix/", model_name, ".eps")
  ggsave(filename, device = "eps", plot = p, width = 8, height = 6, dpi = 300)
}

for (model_name in names(sunab_list)) {
  model <- get(model_name)  # retrieve the model object
  tidy_model <- broom::tidy(model)  # extract coefficients
  
  event_df <- tidy_model %>%
    filter(grepl("TTTyear2::", term)) %>%
      mutate(
        time = as.numeric(gsub("TTTyear2::", "", term)),
      ci_lower = estimate - 1.96 * std.error,
      ci_upper = estimate + 1.96 * std.error
    )
  event_df <- bind_rows(
    event_df,
    tibble(time = -1, estimate = 0, ci_lower = 0, ci_upper = 0)
  )
  
  # Identify group and pull settings
  current_group <- names(Filter(function(x) model_name %in% x, model_groups))[1]
  y_settings <- group_y_axis[[current_group]]
  
  # Plot
  p <- ggplot(event_df, aes(x = time, y = estimate)) +
    geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper), fill = "gray60") +
    geom_line(color = "black") +
    geom_point(color = "black") +
    geom_vline(xintercept = -1, linetype = "dashed", color = "black") +
    geom_hline(yintercept = 0, linetype = "solid", color = "black") + 
    labs(
      title = sunab_list[[model_name]],
      x = "Event Time",
      y = "Estimate and 95% CI"
    ) +
    My_Theme+ 
    scale_y_continuous(breaks = y_settings$breaks,
                       labels = scales::label_number(accuracy = 0.01)) +
    coord_cartesian(ylim = y_settings$limits)  +
    theme(plot.title = element_text(hjust = 0.5))
  
  
  filename <- paste0("Supplementary Appendix/", model_name, ".eps")
  ggsave(filename, device = "eps", plot = p, width = 8, height = 6, dpi = 300)
}

for (model_name in names(closed_list)) {
  model <- get(model_name)  # retrieve the model object
  tidy_model <- broom::tidy(model)  # extract coefficients
  
  event_df <- tidy_model %>%
    filter(grepl("TTT.closedyear2::", term)) %>%
    mutate(
      time = as.numeric(gsub("TTT.closedyear2::", "", term)),
      ci_lower = estimate - 1.96 * std.error,
      ci_upper = estimate + 1.96 * std.error
    )
  event_df <- bind_rows(
    event_df,
    tibble(time = -1, estimate = 0, ci_lower = 0, ci_upper = 0)
  )
  
  # Identify group and pull settings
  current_group <- names(Filter(function(x) model_name %in% x, model_groups))[1]
  y_settings <- group_y_axis[[current_group]]
  
  # Plot
  p <- ggplot(event_df, aes(x = time, y = estimate)) +
    geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper), fill = "gray60") +
    geom_line(color = "black") +
    geom_point(color = "black") +
    geom_vline(xintercept = -1, linetype = "dashed", color = "black") +
    geom_hline(yintercept = 0, linetype = "solid", color = "black") + 
    labs(
      title = closed_list[[model_name]],
      x = "Event Time",
      y = "Estimate and 95% CI"
    ) +
    My_Theme+ 
    scale_y_continuous(breaks = y_settings$breaks,
                       labels = scales::label_number(accuracy = 0.01)) +
    coord_cartesian(ylim = y_settings$limits)  +
    theme(plot.title = element_text(hjust = 0.5))
  
  
  filename <- paste0("Supplementary Appendix/", model_name, ".eps")
  ggsave(filename, device = "eps", plot = p, width = 8, height = 6, dpi = 300)
}


################################################################
#Event studies with varying excluded years 
count_control_obs <- function(model, data, ref_years, treat_var = "first.treat.year", unit_var = "huc12") {
  # Get estimation sample
  obs_model <- obs(model)
  
  # Add row ID if needed
  if (!"row_id" %in% names(data)) {
    data$row_id <- seq_len(nrow(data))
  }
  
  # Filter data to estimation sample
  filtered_data <- data[data$row_id %in% obs_model, ]
  
  # Subset to control units (i.e., units whose treatment year is in ref_years)
  control_units <- filtered_data[filtered_data[[treat_var]] %in% ref_years, ]
  
  # Return number of unique HUC12s in control group
  return(nrow(control_units))
}

ref_list <- list(
  c(9999),                      # aes1
  c(2018),                      # aes2
  2016:2018,                    # aes3
  2014:2018,                    # aes4
  2009:2018                     # aes5
)


#run models with varying cohorts 
#ammonia all huc
aes1<-fixest::feols(ammonia_filtered~sunab(first.treat.year, TTTyear2, att=TRUE) + ..ctrl_n | ..huc8_fe, 
                    data=merge, cluster="huc8year")


#ammonia ever treated. last year as control. 
aes2<-fixest::feols(ammonia_filtered~sunab(first.treat.year, TTTyear2, ref.c=2018, att=TRUE) + ..ctrl_n | ..huc8_fe, 
                    data=merge, cluster="huc8year", subset=merge$EverTreated==1 & merge$Year<2018)

#last 3 years as control. 
aes3<-fixest::feols(ammonia_filtered~sunab(first.treat.year, TTTyear2, ref.c=c(2016, 2017, 2018), att=TRUE) + ..ctrl_n | ..huc8_fe, 
                    data=merge, cluster="huc8year", subset=merge$EverTreated==1 & merge$Year<2016)


#last 5 years as control
aes4<-fixest::feols(ammonia_filtered~sunab(first.treat.year, TTTyear2, ref.c=c(2014, 2015, 2016, 2017, 2018), att=TRUE) + ..ctrl_n | ..huc8_fe, 
                    data=merge, cluster="huc8year", subset=merge$EverTreated==1 & merge$Year<2014)

#last 10 years as control
aes5<-fixest::feols(ammonia_filtered~sunab(first.treat.year, TTTyear2, ref.c=c(2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018), att=TRUE) + ..ctrl_n | ..huc8_fe, 
                    data=merge, cluster="huc8year", subset=merge$EverTreated==1 & merge$Year<2009)
# TKN all HUCs
tes1 <- fixest::feols(TKN_filtered ~ sunab(first.treat.year, TTTyear2, att=TRUE) + ..ctrl_n | ..huc8_fe, 
                      data=merge, cluster="huc8year")

# TKN ever treated, last year as control (2018)
tes2 <- fixest::feols(TKN_filtered ~ sunab(first.treat.year, TTTyear2, ref.c=2018, att=TRUE) + ..ctrl_n | ..huc8_fe, 
                      data=merge, cluster="huc8year", subset=merge$EverTreated == 1 & merge$Year<2018)


# TKN last 3 years as control
tes3 <- fixest::feols(TKN_filtered ~ sunab(first.treat.year, TTTyear2, ref.c=c(2016,2017,2018), att=TRUE) + ..ctrl_n | ..huc8_fe, 
                      data=merge, cluster="huc8year", subset=merge$EverTreated == 1 & merge$Year<2016)

# TKN last 5 years as control
tes4 <- fixest::feols(TKN_filtered ~ sunab(first.treat.year, TTTyear2, ref.c=c(2014, 2015, 2016, 2017, 2018), att=TRUE) + ..ctrl_n | ..huc8_fe, 
                      data=merge, cluster="huc8year", subset=merge$EverTreated == 1 & merge$Year<2014)

# TKN last 10 years as control
tes5 <- fixest::feols(TKN_filtered ~ sunab(first.treat.year, TTTyear2, ref.c=c(2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018), att=TRUE) + ..ctrl_n | ..huc8_fe, 
                      data=merge, cluster="huc8year", subset=merge$EverTreated == 1 & merge$Year<2009)

# Phosphorus all HUCs
pes1 <- fixest::feols(TP_filtered ~ sunab(first.treat.year, TTTyear2, att=TRUE) + ..ctrl_p | ..huc8_fe, 
                      data=merge, cluster="huc8year")

# Phosphorus ever treated, last year as control (2018)
pes2 <- fixest::feols(TP_filtered ~ sunab(first.treat.year, TTTyear2, ref.c=2018, att=TRUE) + ..ctrl_p | ..huc8_fe, 
                      data=merge, cluster="huc8year", subset=merge$EverTreated == 1 & merge$Year<2018)

# Phosphorus last 3 years as control
pes3 <- fixest::feols(TP_filtered ~ sunab(first.treat.year, TTTyear2, ref.c=c(2016,2017,2018), att=TRUE) + ..ctrl_p | ..huc8_fe, 
                      data=merge, cluster="huc8year", subset=merge$EverTreated == 1 & merge$Year<2016)

# Phosphorus last 5 years as control
pes4 <- fixest::feols(TP_filtered ~ sunab(first.treat.year, TTTyear2, ref.c=c(2014, 2015, 2016, 2017, 2018), att=TRUE) + ..ctrl_p | ..huc8_fe, 
                      data=merge, cluster="huc8year", subset=merge$EverTreated == 1 & merge$Year<2014)

# Phosphorus last 10 years as control
pes5 <- fixest::feols(TP_filtered ~ sunab(first.treat.year, TTTyear2, ref.c=c(2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018), att=TRUE) + ..ctrl_p | ..huc8_fe, 
                      data=merge, cluster="huc8year", subset=merge$EverTreated == 1 & merge$Year<2009)


#output ATT tables 

#table ammonia 
etable(aes1, aes2, aes3, aes4, aes5,
       keep=c("ATT"),
       drop.section = "fixef",
       fitstat=c("my", "n", "ar2"),
       style.tex=fixest::style.tex("aer"),
       tex=TRUE,
       file="Supplementary Appendix/ammonia_sa_robust.tex",
       replace=TRUE,
       fixef_sizes=TRUE, 
       fixef_sizes.simplify = TRUE,
       extralines=list("__Control Group"=c("NT", "NYT 2018", "NYT 2016-2018", "NYT 2014-2018", "NYT 2009-2018"),
                       "__# HUC12" = sapply(list(aes1, aes2, aes3, aes4, aes5), count_huc12, data = merge),
                       "__# Control HUC12" = sapply(seq_along(ref_list), function(i) {
                         count_control_huc12(model = list(aes1, aes2, aes3, aes4, aes5)[[i]],
                                             data = merge,
                                             ref_years = ref_list[[i]])
                       }),
                       "__# Control Obs" = sapply(seq_along(ref_list), function(i) {
                         count_control_obs(model = list(aes1, aes2, aes3, aes4, aes5)[[i]],
                                           data=merge,
                                           ref_years = ref_list[[i]])
                       }))
)
#table TKN 

etable(tes1, tes2, tes3, tes4, tes5,
       keep=c("ATT"),
       drop.section = "fixef",
       fitstat=c("my", "n", "ar2"),
       style.tex=fixest::style.tex("aer"),
       tex=TRUE,
       file="Supplementary Appendix/tkn_sa_robust.tex",
       replace=TRUE,
       fixef_sizes=TRUE, 
       fixef_sizes.simplify = TRUE,
       extralines=list("__Control Group"=c("NT", "NYT 2018", "NYT 2016-2018", "NYT 2014-2018", "NYT 2009-2018"),
                       "__# HUC12" = sapply(list(tes1, tes2, tes3, tes4, tes5), count_huc12, data = merge),
                       "__# Control HUC12" = sapply(seq_along(ref_list), function(i) {
                         count_control_huc12(model = list(tes1, tes2, tes3, tes4, tes5)[[i]],
                                             data = merge,
                                             ref_years = ref_list[[i]])
                       }),
                       "__# Control Obs" = sapply(seq_along(ref_list), function(i) {
                         count_control_obs(model = list(tes1, tes2, tes3, tes4, tes5)[[i]],
                                           data=merge,
                                           ref_years = ref_list[[i]])
                       }))
)


#Table for phosphorus
etable(pes1, pes2, pes3, pes4, pes5,
       keep=c("ATT"),
       drop.section = "fixef",
       fitstat=c("my", "n", "ar2"),
       style.tex=fixest::style.tex("aer"),
       tex=TRUE,
       file="Supplementary Appendix/phosphorus_sa_robust.tex",
       replace=TRUE,
       fixef_sizes=TRUE, 
       fixef_sizes.simplify = TRUE,
       extralines=list("__Control Group"=c("NT", "NYT 2018", "NYT 2016-2018", "NYT 2014-2018", "NYT 2009-2018"),
                       "__# HUC12" = sapply(list(pes1, pes2, pes3, pes4, pes5), count_huc12, data = merge),
                       "__# Control HUC12" = sapply(seq_along(ref_list), function(i) {
                         count_control_huc12(model = list(pes1, pes2, pes3, pes4, pes5)[[i]],
                                             data = merge,
                                             ref_years = ref_list[[i]])
                       }),
                       "__# Control Obs" = sapply(seq_along(ref_list), function(i) {
                         count_control_obs(model = list(pes1, pes2, pes3, pes4, pes5)[[i]],
                                           data=merge,
                                           ref_years = ref_list[[i]])
                       }))
)



###################################################################################
#Heterogeneity analysis 


#ag categories 
ever_treated <- subset(merge, EverTreated == 1)
quantile(ever_treated$cropland, probs = c(.33, .67))
ever_treated$cropland_quartile <- cut(ever_treated$cropland, breaks = 3)
ag_q1 <- subset(ever_treated, cropland <= 0.25)
ag_q2 <- subset(ever_treated, cropland > 0.25 & cropland <= 0.50 )
ag_q3 <- subset(ever_treated, cropland > 0.50 & cropland <=0.75 )
ag_q4 <- subset(ever_treated, cropland > 0.75 )

tkn_q2 <- subset(ag_q2, !is.na(TKN_filtered))
tkn_q4 <- subset(ag_q4, !is.na(TKN_filtered))
length(unique(tkn_q2$huc12))
length(unique(tkn_q4$huc12))


#manually add the means in :(
summary(ag_q1$ammonia_filtered)
summary(ag_q2$ammonia_filtered)
summary(ag_q3$ammonia_filtered)
summary(ag_q4$ammonia_filtered)

summary(ag_q1$TP_filtered)
summary(ag_q2$TP_filtered)
summary(ag_q3$TP_filtered)
summary(ag_q4$TP_filtered)

summary(ag_q1$TKN_filtered)
summary(ag_q2$TKN_filtered)
summary(ag_q3$TKN_filtered)
summary(ag_q4$TKN_filtered)

#Dependent variable mean (ammonia) & 0.065 & 0.065 & .065 & 0.096 & 0.096 & 0.096 & 0.122 & 0.122 & 0.122 & 0.250 & 0.250 & 0.250 \\ 
#Dependent variable mean (phosphorus) & 0.094 & 0.094 & 0.094 & 0.095 & 0.095& 0.095  & 0.153  & 0.153 & 0.153 & 0.232   & 0.232  & 0.232 \\ 
#Dependent variable mean (TKN) & 0.310 & 0.310 & 0.310 & 0.370 & 0.370 & 0.370 & 0.571 & 0.571 & 0.571 & 0.710 & 0.710 & 0.710 \\ 
for (q in 1:4){
  data_q <- get(paste0("ag_q", q))
  model <- feols(fml=ammonia_filtered~ihs_restored +..ctrl_n| ..huc8_fe,
               cluster="huc8year", data=data_q)
  assign(paste0("a1_q", q), model)
  model <- feols(fml=ammonia_filtered~ihs_restored+..ctrl_n| ..huc4year_fe,
               cluster="huc8year", data=data_q)
  assign(paste0("a2_q", q), model)
  model <- feols(fml=ammonia_filtered~ihs_restored+..ctrl_n| ..huc12_fe,
                cluster="huc8year", data=data_q)
  assign(paste0("a3_q", q), model)
  model <- feols(fml=TP_filtered~ihs_restored +..ctrl_p| ..huc8_fe,
        cluster="huc8year", data=data_q)
  assign(paste0("p1_q", q), model)
  model <- feols(fml=TP_filtered~ihs_restored+..ctrl_p| ..huc4year_fe,
        cluster="huc8year", data=data_q)
  assign(paste0("p2_q", q), model)
  model <- feols(fml=TP_filtered~ihs_restored+..ctrl_p| ..huc12_fe,
        cluster="huc8year", data=data_q)
  assign(paste0("p3_q", q), model)
  model <- feols(fml=TKN_filtered~ihs_restored +..ctrl_n| ..huc8_fe,
        cluster="huc8year", data=data_q)
  assign(paste0("n1_q", q), model)
  model <- feols(fml=TKN_filtered~ihs_restored+..ctrl_n| ..huc4year_fe,
        cluster="huc8year", data=data_q)
  assign(paste0("n2_q", q), model)
  model <- feols(fml=TKN_filtered~ihs_restored+..ctrl_n| ..huc12_fe,
        cluster="huc8year", data=data_q)
  assign(paste0("n3_q", q), model)
  
}

models_a <- list(a1_q1, a2_q1, a3_q1, a1_q2, a2_q2, a3_q2, a1_q3, a2_q3, a3_q3, a1_q4, a2_q4, a3_q4)
models_p <- list(p1_q1, p2_q1, p3_q1, p1_q2, p2_q2, p3_q2, p1_q3, p2_q3, p3_q3, p1_q4, p2_q4, p3_q4)
models_n <- list(n1_q1, n2_q1, n3_q1, n1_q2, n2_q2, n3_q2, n1_q3, n2_q3, n3_q3, n1_q4, n2_q4, n3_q4)

fixest::etable(models_a,
               depvar = FALSE, drop.section="fixef",
               tex=TRUE, digits=3, fitstat=c("n", "ar2"), 
               keep=c("IHS"),
               file="Supplementary Appendix/ammonia_ag_percentile.tex",
               replace = TRUE,               
               postprocess.tex = change_column_widths,
               style.tex=fixest::style.tex("aer", line.bottom = "",
                                           var.title = "$\\textbf{Panel A: ammonia}$"),
               extralines=list("__ "=c(" ", " ", " ", " ", " ", " ", "" , "" , "", "" , "" , ""),
                               "^^  "=c(" ", " ", " ", " ", " ", " ", "" , "" , "", "" , "" , "")))

fixest::etable(n1_q1, n2_q1, n3_q1, n1_q2, n2_q2, n3_q2, n1_q3, n2_q3, n3_q3, n1_q4, n2_q4, n3_q4,
               depvar = FALSE, drop.section="fixef",
               tex=TRUE, digits=3, fitstat=c("n", "ar2"), 
               keep=c("IHS"),
               file="Supplementary Appendix/tkn_ag_percentile.tex",
               replace = TRUE,               
               postprocess.tex = both,
               style.tex=fixest::style.tex("aer", line.top = "", line.bottom = "",
                                           var.title = "$\\textbf{Panel B: TKN}$"),
               extralines=list("__ "=c(" ", " ", " ", " ", " ", " ", "" , "" , "", "" , "" , ""),
                               "^^  "=c(" ", " ", " ", " ", " ", " ", "" , "" , "", "" , "" , "")))


fixest::etable(p1_q1, p2_q1, p3_q1, p1_q2, p2_q2, p3_q2, p1_q3, p2_q3, p3_q3, p1_q4, p2_q4, p3_q4,
               depvar = FALSE, drop.section = "fixef",
               tex=TRUE, digits=3, fitstat=c("n", "ar2"), 
               keep=c("IHS"),
               file="Supplementary Appendix/phosphorous_ag_percentile.tex",
               replace = TRUE,
               style.tex=fixest::style.tex("aer", line.top ="",
                                           fixef.title = "\\midrule",
                                           var.title = "$\\textbf{Panel C: phosphorus}$"),
               postprocess.tex = both,
               extralines=list("__"=c(" ", " ", " ", " "," ", " ", "" ,"" , "", " ", "" ,"" ),
                               "__Ever treated"=c("$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$","$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$"),
                               "__Controls"=c("$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$","$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$"),
                               "__"=c(" ", " ", " ", " "," ", " ", " "," ", " ", " "," ", " "),
                               "__HUC8 fixed effects"=c("$\\checkmark$",  "$\\checkmark$"," ", "$\\checkmark$", "$\\checkmark$"," ", "$\\checkmark$", "$\\checkmark$"," ", "$\\checkmark$", "$\\checkmark$"," "),
                               "__HUC12 fixed effects"=c(" "," ", "$\\checkmark$", " ", " ", "$\\checkmark$", " ", " ", "$\\checkmark$", " ", " ", "$\\checkmark$"),
                               "__Year fixed effects"=c("$\\checkmark$"," ", "$\\checkmark$", "$\\checkmark$"," ","$\\checkmark$", "$\\checkmark$"," ","$\\checkmark$", "$\\checkmark$"," ","$\\checkmark$"),
                               "__HUC4 x year fixed effects"=c(" ", "$\\checkmark$", " "," ", "$\\checkmark$"," "," ", "$\\checkmark$"," "," ", "$\\checkmark$"," "),
                               "__Month fixed effects"=c("$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$","$\\checkmark$", "$\\checkmark$", "$\\checkmark$","$\\checkmark$", "$\\checkmark$", "$\\checkmark$","$\\checkmark$", "$\\checkmark$")))



#weather interactions 
#first define our linear hypothesis tests 
test_coefficients <- function(model) {
  # Specify the hypotheses to test
glht_test <-  glht(model, linfct = c("ihs_restored + ihs_restored:PPT_binLow = 0",
                        "ihs_restored + ihs_restored:PPT_binMedium = 0",
                        "ihs_restored + ihs_restored:PPT_binHigh = 0",
                        "ihs_restored + ihs_restored:PPT_binVeryHigh = 0",
                        "ihs_restored + ihs_restored:Temp_binLow = 0",
                        "ihs_restored + ihs_restored:Temp_binMedium = 0",
                        "ihs_restored + ihs_restored:Temp_binHigh = 0",
                        "ihs_restored + ihs_restored:Temp_binVeryHigh = 0") )  
  # Return the summary of the test results
  return(summary(glht_test))
}
a1_w<-feols(fml=ammonia_filtered~ihs_restored+ ..int_n | ..huc8_fe, 
          cluster="huc8year", data=merge, subset=merge$EverTreated==1)
a2_w<-feols(fml=ammonia_filtered~ihs_restored+ ..int_n| ..huc4year_fe,
          cluster="huc8year", data=merge, subset=merge$EverTreated==1)

a3_w<-feols(fml=ammonia_filtered~ihs_restored+..int_n| ..huc12_fe,
          cluster="huc8year", data=merge, subset=merge$EverTreated==1)

#p
p1_w<-feols(fml=TP_filtered~ihs_restored+ ..int_p | ..huc8_fe, 
          cluster="huc8year", data=merge, subset=merge$EverTreated==1)

p2_w<-feols(fml=TP_filtered~ihs_restored+ ..int_p| ..huc4year_fe,
          cluster="huc8year", data=merge, subset=merge$EverTreated==1)

p3_w<-feols(fml=TP_filtered~ihs_restored+..int_p| ..huc12_fe,
          cluster="huc8year", data=merge, subset=merge$EverTreated==1)

test_coefficients(p1_w)
test_coefficients(p2_w)
test_coefficients(p3_w)

#tkn 
n1_w<-feols(fml=TKN_filtered~ihs_restored+ ..int_n | ..huc8_fe, 
          cluster="huc8year", data=merge, subset=merge$EverTreated==1)

n2_w<-feols(fml=TKN_filtered~ihs_restored+ ..int_n| ..huc4year_fe,
          cluster="huc8year", data=merge, subset=merge$EverTreated==1)

n3_w<-feols(fml=TKN_filtered~ihs_restored+..int_n| ..huc12_fe,
          cluster="huc8year", data=merge, subset=merge$EverTreated==1)


fixest::etable(a1_w, a2_w, a3_w, n1_w, n2_w, n3_w, p1_w, p2_w, p3_w,
               depvar = FALSE, drop.section="fixef",
               tex=TRUE, digits=3, fitstat=c("my", "n", "ar2"), 
               keep=c("IHS"),
               file="Supplementary Appendix/ammonia_weather.tex",
               replace = TRUE,               
               postprocess.tex = change_column_widths,
               style.tex=fixest::style.tex("aer", line.bottom = ""),
               extralines=list("__"=c(" ", " ", " ", " ", " ", " ", " ", " ", " "),
                               "__Ever treated"=c("$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$","$\\checkmark$", "$\\checkmark$", "$\\checkmark$"),
                               "__Controls"=c("$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$","$\\checkmark$", "$\\checkmark$", "$\\checkmark$"),
                               "__"=c(" ", " ", " "," ", " ", " "," ", " ", " "),
                               "__HUC8 fixed effects"=c("$\\checkmark$",  "$\\checkmark$"," ","$\\checkmark$",  "$\\checkmark$"," ","$\\checkmark$",  "$\\checkmark$"," "),
                               "__HUC12 fixed effects"=c(" "," ", "$\\checkmark$"," "," ", "$\\checkmark$"," "," ", "$\\checkmark$"),
                               "__Year fixed effects"=c("$\\checkmark$"," ", "$\\checkmark$","$\\checkmark$"," ", "$\\checkmark$","$\\checkmark$"," ", "$\\checkmark$"),
                               "__HUC4 x year fixed effects"=c(" ", "$\\checkmark$", " "," ", "$\\checkmark$", " "," ", "$\\checkmark$", " "),
                               "__Month fixed effects"=c("$\\checkmark$", "$\\checkmark$", "$\\checkmark$","$\\checkmark$", "$\\checkmark$", "$\\checkmark$","$\\checkmark$", "$\\checkmark$", "$\\checkmark$")))
fixest::etable(n1_w, n2_w, n3_w, 
               depvar = FALSE, drop.section="fixef",
               tex=TRUE, digits=3, fitstat=c("my", "n", "ar2"), 
               keep=c("IHS"),
               file="Supplementary Appendix/tkn_weather.tex",
               replace = TRUE,               
               postprocess.tex = both,
               style.tex=fixest::style.tex("aer", line.top = "", line.bottom = "",
                                           var.title = "$\\textbf{Panel B: TKN}$"),
               extralines=list("__ "=c(" ", " ", " "),
                               "^^  "=c(" ", " ", " ")))


fixest::etable(p1_w, p2_w, p3_w,
               depvar = FALSE, drop.section = "fixef",
               tex=TRUE, digits=3, fitstat=c("my", "n", "ar2"), 
               keep=c("IHS"),
               file="Supplementary Appendix/phosphorus_weather.tex",
               replace = TRUE,
               style.tex=fixest::style.tex("aer", line.top ="",
                                           fixef.title = "\\midrule",
                                           var.title = "$\\textbf{Panel C: phosphorus}$"),
               postprocess.tex = both,
               extralines=list("__"=c(" ", " ", " "),
                               "__Ever treated"=c("$\\checkmark$", "$\\checkmark$", "$\\checkmark$"),
                               "__Controls"=c("$\\checkmark$", "$\\checkmark$", "$\\checkmark$"),
                               "__"=c(" ", " ", " "),
                               "__HUC8 fixed effects"=c("$\\checkmark$",  "$\\checkmark$"," "),
                               "__HUC12 fixed effects"=c(" "," ", "$\\checkmark$"),
                               "__Year fixed effects"=c("$\\checkmark$"," ", "$\\checkmark$"),
                               "__HUC4 x year fixed effects"=c(" ", "$\\checkmark$", " "),
                               "__Month fixed effects"=c("$\\checkmark$", "$\\checkmark$", "$\\checkmark$")))

###################################################################################
#Alternative nutrient outcomes 

#need to reload the data so we have all the observations of the alternative nutrients 
load("Data/Main/main_huc12_panel.RData")


#keep only relevant years 
merge<-subset(merge, as.numeric(Year)>=1985)


#Ammonia unfiltered
a1<-feols(fml=ammonia_unfiltered~ihs_restored+..ctrl_n| ..huc8_fe, 
          cluster="huc8year", data=merge)

a2<-feols(fml=ammonia_unfiltered~ihs_restored+ihs_upstream_restored+..ctrl_n| ..huc8_fe,
          cluster="huc8year", data=merge)

a3<-feols(fml=ammonia_unfiltered~ihs_restored+..ctrl_n| ..huc8_fe, 
          cluster="huc8year", data=merge, subset=merge$EverTreated==1)

a4<-feols(fml=ammonia_unfiltered~ihs_restored+ihs_upstream_restored+..ctrl_n| ..huc8_fe,
          cluster="huc8year", data=merge, subset=merge$EverTreated==1)

#TKN unfiltered
tkn1<-feols(fml=ammonia_unfiltered~ihs_restored+..ctrl_n| ..huc8_fe, 
            cluster="huc8year", data=merge)

tkn2<-feols(fml=ammonia_unfiltered~ihs_restored+ihs_upstream_restored+..ctrl_n| ..huc8_fe,
            cluster="huc8year", data=merge)

tkn3<-feols(fml=ammonia_unfiltered~ihs_restored+..ctrl_n| ..huc8_fe, 
            cluster="huc8year", data=merge, subset=merge$EverTreated==1)

tkn4<-feols(fml=ammonia_unfiltered~ihs_restored+ihs_upstream_restored+..ctrl_n| ..huc8_fe,
            cluster="huc8year", data=merge, subset=merge$EverTreated==1)


#nitrogen unfiltered
iu1<-feols(fml=TN_unfiltered~ihs_restored+..ctrl_n| ..huc8_fe, 
           cluster="huc8year", data=merge)

iu2<-feols(fml=TN_unfiltered~ihs_restored+ihs_upstream_restored+..ctrl_n| ..huc8_fe,
           cluster="huc8year", data=merge)

iu3<-feols(fml=TN_unfiltered~ihs_restored+..ctrl_n| ..huc8_fe, 
           cluster="huc8year", data=merge, subset=merge$EverTreated==1)

iu4<-feols(fml=TN_unfiltered~ihs_restored+ihs_upstream_restored+..ctrl_n| ..huc8_fe,
           cluster="huc8year", data=merge, subset=merge$EverTreated==1)



fixest::etable(a1, a2, a3, a4,
               depvar = FALSE, drop.section = "fixef",
               tex=TRUE, digits=3, fitstat=c("my", "n", "ar2"), 
               file = "Supplementary Appendix/ammonia_unfiltered_twfe.tex",
               keep=c("IHS"),  powerBelow=-9,
               replace = TRUE,
               postprocess.tex = change_column_widths,
               style.tex=fixest::style.tex("aer", line.bottom = "", var.title = "$\\textbf{Ammonia (unfiltered)}$"))

fixest::etable(tkn1, tkn2,tkn3, tkn4,
               depvar = FALSE, drop.section = "fixef",
               tex=TRUE, digits=3, fitstat=c("my", "n", "ar2"), 
               file = "Supplementary Appendix/tkn_unfiltered_twfe.tex",
               keep=c("IHS"),  powerBelow=-9,
               replace = TRUE,
               postprocess.tex = change_column_widths,
               style.tex=fixest::style.tex("aer", line.bottom = "", var.title = "$\\textbf{TKN (unfiltered)}$"))

fixest::etable(iu1, iu2, iu3, iu4, 
               depvar = FALSE, drop.section = "fixef",
               tex=TRUE, digits=3, fitstat=c("my", "n", "ar2"),
               file = "Supplementary Appendix/nitrogen_unfiltered_twfe.tex",
               keep=c("IHS"),  powerBelow=-9,
               replace = TRUE,
               postprocess.tex = change_column_widths,
               style.tex=fixest::style.tex("aer", line.bottom = "", var.title = "$\\textbf{Nitrogen (unfiltered)}$"))

#Phosphorus nutrients
#phosphorus unfiltered
t1<-feols(fml=TP_unfiltered~ihs_restored+..ctrl_p| ..huc8_fe, 
          cluster="huc8year", data=merge)

t2<-feols(fml=TP_unfiltered~ihs_restored+ihs_upstream_restored+..ctrl_p| ..huc8_fe,
          cluster="huc8year", data=merge)

t3<-feols(fml=TP_unfiltered~ihs_restored+..ctrl_p| ..huc8_fe, 
          cluster="huc8year", data=merge, subset=merge$EverTreated==1)

t4<-feols(fml=TP_unfiltered~ihs_restored+ihs_upstream_restored+..ctrl_p| ..huc8_fe,
          cluster="huc8year", data=merge, subset=merge$EverTreated==1)


#orthophosphate filtered
o1<-feols(fml=orthophosphate_filtered~ihs_restored+..ctrl_p| ..huc8_fe, 
          cluster="huc8year", data=merge)

o2<-feols(fml=orthophosphate_filtered~ihs_restored+ihs_upstream_restored+..ctrl_p| ..huc8_fe,
          cluster="huc8year", data=merge)

o3<-feols(fml=orthophosphate_filtered~ihs_restored+..ctrl_p| ..huc8_fe, 
          cluster="huc8year", data=merge, subset=merge$EverTreated==1)

o4<-feols(fml=orthophosphate_filtered~ihs_restored+ihs_upstream_restored+..ctrl_p| ..huc8_fe,
          cluster="huc8year", data=merge, subset=merge$EverTreated==1)

#orthophosphate unfiltered
ou1<-feols(fml=orthophosphate_unfiltered~ihs_restored+..ctrl_p| ..huc8_fe, 
           cluster="huc8year", data=merge)

ou2<-feols(fml=orthophosphate_unfiltered~ihs_restored+ihs_upstream_restored+..ctrl_p| ..huc8_fe,
           cluster="huc8year", data=merge)

ou3<-feols(fml=orthophosphate_unfiltered~ihs_restored+..ctrl_p| ..huc8_fe, 
           cluster="huc8year", data=merge, subset=merge$EverTreated==1)

ou4<-feols(fml=orthophosphate_unfiltered~ihs_restored+ihs_upstream_restored+..ctrl_p| ..huc8_fe,
           cluster="huc8year", data=merge, subset=merge$EverTreated==1)

fixest::etable(t1, t2, t3, t4,
               depvar = FALSE, drop.section = "fixef",
               tex=TRUE, digits=3, fitstat=c("my", "n", "ar2"), 
               file = "Supplementary Appendix/phosphorus_unfiltered_twfe.tex",
               keep=c("IHS"), powerBelow=-9,
               replace = TRUE,
               postprocess.tex = change_column_widths,
               style.tex=fixest::style.tex("aer", line.bottom = "", var.title = "$\\textbf{Phosphorus (unfiltered)}$"))

fixest::etable(o1, o2, o3, o4,
               depvar = FALSE, drop.section = "fixef",
               tex=TRUE, digits=3, fitstat=c("my", "n", "ar2"), 
               file = "Supplementary Appendix/orthophosphate_filtered_twfe.tex",
               keep=c("IHS"), powerBelow=-9,
               replace = TRUE,
               postprocess.tex = change_column_widths,
               style.tex=fixest::style.tex("aer", line.bottom = "", var.title = "$\\textbf{Orthophosphate (filtered)}$"))


fixest::etable(ou1, ou2, ou3, ou4,
               depvar = FALSE, drop.section = "fixef",
               tex=TRUE, digits=3, fitstat=c("my", "n", "ar2"), 
               file = "Supplementary Appendix/orthophosphate_unfiltered_twfe.tex",
               keep=c("IHS"), powerBelow=-9,
               replace = TRUE,
               postprocess.tex = change_column_widths,
               style.tex=fixest::style.tex("aer", line.bottom = "", var.title = "$\\textbf{Orthophosphate (unfiltered)}$"))


#Event studies 
#Ammonia unfiltered
aues1<-fixest::feols(ammonia_unfiltered~sunab(first.treat.year, TTTyear2) + ..ctrl_n | ..huc8_fe, 
                     data=merge, cluster="huc8year", subset=merge$EverTreated==1)

aues2<-fixest::feols(ammonia_unfiltered~sunab(first.treat.year, TTTyear2) + ..ctrl_n  | ..huc8_fe, 
                     data=merge, cluster="huc8year")

#TKN unfiltered
tknues1<-fixest::feols(TKN_unfiltered~sunab(first.treat.year, TTTyear2) + ..ctrl_n  | ..huc8_fe, 
                       data=merge, cluster="huc8year", subset=merge$EverTreated==1)

tknues2<-fixest::feols(TKN_unfiltered~sunab(first.treat.year, TTTyear2) + ..ctrl_n  | ..huc8_fe, 
                       data=merge, cluster="huc8year")

#TN unfiltered
tnues1<-fixest::feols(TN_unfiltered~sunab(first.treat.year, TTTyear2) + ..ctrl_n  | ..huc8_fe, 
                      data=merge, cluster="huc8year", subset=merge$EverTreated==1)

tnues2<-fixest::feols(TN_unfiltered~sunab(first.treat.year, TTTyear2) + ..ctrl_n  | ..huc8_fe, 
                      data=merge, cluster="huc8year")

#TP unfiltered 
tpu1<-fixest::feols(TP_unfiltered~sunab(first.treat.year, TTTyear2) + ..ctrl_p | ..huc8_fe, 
                    data=merge, cluster="huc8year", subset=merge$EverTreated==1)
tpu2<-fixest::feols(TP_unfiltered~sunab(first.treat.year, TTTyear2) + ..ctrl_p | ..huc8_fe, 
                    data=merge, cluster="huc8year")

#orthophosphate_filtered
o1<-fixest::feols(orthophosphate_filtered~sunab(first.treat.year, TTTyear2) + ..ctrl_p | ..huc8_fe, 
                  data=merge, cluster="huc8year", subset=merge$EverTreated==1)
o2<-fixest::feols(orthophosphate_filtered~sunab(first.treat.year, TTTyear2) + ..ctrl_p | ..huc8_fe, 
                  data=merge, cluster="huc8year")

#orthophosphate_unfiltered
ou1<-fixest::feols(orthophosphate_unfiltered~sunab(first.treat.year, TTTyear2) + ..ctrl_p | ..huc8_fe, 
                   data=merge, cluster="huc8year", subset=merge$EverTreated==1)
ou2<-fixest::feols(orthophosphate_unfiltered~sunab(first.treat.year, TTTyear2) + ..ctrl_p | ..huc8_fe, 
                   data=merge, cluster="huc8year")


#Plot them 
My_Theme = theme(title=element_text(size=12),
                 text = element_text(color = "black"),
                 strip.text = element_text(color = "black", size=12),
                 strip.background=element_rect(colour=NA,
                                               fill=NA),
                 axis.title.x = element_text(size = 12),
                 axis.text.x = element_text(size = 12),
                 axis.text.y = element_text(size = 12),
                 axis.title.y = element_text(size = 12),
                 legend.position = "none")

  
  
  ggiplot(list("Ammonia (unfiltered) \n Ever-treated HUC12"=aues1, 
                                 "Ammonia (unfiltered) \n All HUC12"= aues2), #ammonia unfiltered
                            multi_style="facet",
                            main=" ",
                            facet_args=list(ncol=2, scales="fixed"),
                            ylab="Estimate and 95% Conf. Int.",
                            pt.join=TRUE,
                            geom_style = 'ribbon',
                            pt.col=c("black", "black","black", "black"), col=c("black","black","black", "black"),
                            xlab="Time to treatment (year)",)+My_Theme
robust1<-ammonia_unfiltered+scale_y_continuous(limits=c(-.1, 0.1), breaks=seq(-0.1, 0.1, by=0.05))
robust1

tkn_unfiltered<-ggiplot(list("TKN (unfiltered) \n Ever-treated HUC12"=tknues1, 
                             "TKN (unfiltered) \n All HUC12"=tknues2), #TKN unfiltered
                        multi_style="facet",
                        facet_args=list(ncol=2),
                        main=" ",
                        ylab="Estimate and 95% Conf. Int.",
                        pt.join=TRUE,
                        geom_style = 'ribbon',
                        pt.col=c("black", "black","black", "black"), col=c("black","black","black", "black"),
                        xlab="Time to treatment (year)",)+My_Theme
robust2<-tkn_unfiltered+scale_y_continuous(limits=c(-.15, 0.15), breaks=seq(-0.15, 0.15, by=0.05),
                                           labels = scales::number_format(accuracy = 0.01))
nitrogen_unfiltered<-ggiplot(list("Nitrogen (unfiltered) \n Ever-treated HUC12" = tnues1,
                                  "Nitrogen (unfiltered) \n All HUC12" =tnues2), #nitrogen unfiltered
                             multi_style="facet",
                             facet_args=list(ncol=2),
                             main=" ",
                             ylab="Estimate and 95% Conf. Int.",
                             pt.join=TRUE,
                             geom_style = 'ribbon',
                             pt.col=c("black", "black","black", "black"), col=c("black","black","black", "black"),
                             xlab="Time to treatment (year)",)+My_Theme
robust3<-nitrogen_unfiltered+scale_y_continuous(limits=c(-1, 1), breaks=seq(-1,1, by=0.5),
                                                labels = scales::number_format(accuracy = 0.1))

#bind them together
combine1<-robust1/robust2/robust3

ggsave(combine1, device = grDevices::cairo_ps, file="Supplementary Appendix/othernutrients_es1.eps",
       width=9, height=12, units="in")

TP_unfiltered<-ggiplot(list("TP (unfiltered) \n Ever-treated HUC12"=tpu1, 
                            "TP (unfiltered) \n All HUC12"= tpu2), #TP unfiltered
                       multi_style="facet",
                       main=" ",
                       facet_args=list(ncol=2, scales="fixed"),
                       ylab="Estimate and 95% Conf. Int.",
                       pt.join=TRUE,
                       geom_style = 'ribbon',
                       pt.col=c("black", "black","black", "black"), col=c("black","black","black", "black"),
                       xlab="Time to treatment (year)",)+My_Theme

TP_unfiltered1<-TP_unfiltered+scale_y_continuous(limits=c(-.1, 0.1), breaks=seq(-0.1, 0.1, by=0.05))

#orthophosphate filtered
o_filtered<-ggiplot(list("Orthophosphate (filtered) \n Ever-treated HUC12"=o1, 
                         "Orthophosphate (filtered) \n All HUC12"= o2), #orthophosphate filtered
                    multi_style="facet",
                    main=" ",
                    facet_args=list(ncol=2, scales="fixed"),
                    ylab="Estimate and 95% Conf. Int.",
                    pt.join=TRUE,
                    geom_style = 'ribbon',
                    pt.col=c("black", "black","black", "black"), col=c("black","black","black", "black"),
                    xlab="Time to treatment (year)",)+My_Theme

o_filtered1<-o_filtered+scale_y_continuous(limits=c(-.1, 0.1), breaks=seq(-0.1, 0.1, by=0.05))

#orthohosphate unfiltered
o_unfiltered<-ggiplot(list("Orthophosphate (unfiltered) \n Ever-treated HUC12"=ou1, 
                           "Orthophosphate (unfiltered) \n All HUC12"= ou2), #orthophosphate filtered
                      multi_style="facet",
                      main=" ",
                      facet_args=list(ncol=2, scales="fixed"),
                      ylab="Estimate and 95% Conf. Int.",
                      pt.join=TRUE,
                      geom_style = 'ribbon',
                      pt.col=c("black", "black","black", "black"), col=c("black","black","black", "black"),
                      xlab="Time to treatment (year)",)+My_Theme
o_unfiltered1<-o_unfiltered+scale_y_continuous(limits=c(-.1, 0.1), breaks=seq(-0.1, 0.1, by=0.05))

combinep<-TP_unfiltered1/o_filtered1/o_unfiltered1
ggsave(combinep, device = grDevices::cairo_ps, file="Supplementary Appendix/p_nutrients_es_sa.eps",
       width=9, height=12, units="in")

